Quantum Monte Carlo Simulation of the two-dimensional ionic Hubbard model 
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(Dated: February 1, 2010) 

The Quantum Monte Carlo simulations of the ionic Hubbard model on a two-dimensional square 
lattice at half filling were performed. The method based on the direct-space, proposed by Suzuki and 
al., Hirsch and al., was used. Cycles of increasing and decreasing values of the Coulomb interaction 
U were performed for fixed temperature (kT = 0.01). Results indicate that, at low temperature, 
the two insulator phases are separated by a metallic phase for weak to intermediate values of the 
staggered potential A. For large Coulomb repulsion the system is in a Mott insulator with an 
antiferromagnetism order. On increasing and decreasing the Coulomb interaction U the metal-Mott 
insulator transition shows an hysteresis phenomenon while the metal-band insulator transition is 
continue. For large A it seems that the metallic region shrinks to a single metallic point. However, 
the band insulator to the Mott insulator transition is not direct for the studied model. A phase 
diagram is drawn for the temperature kT = 0.01. For A = 0.5 cycles of increasing and decreasing 
temperature were programmed for different values of the Coulomb interaction U . A behaviour 
change appears for U ~ 1.75. This suggests that a crossover line divides the metallic region of the 
phase diagram. 

PACS numbers: 71.27.+a, 71.10.Fd, 71.30.+h 



I. INTRODUCTION 

Recently some theoretical and numerical studies were 
published which investigate the metal-insulator transi- 
tions and the transition between the two insulator phases 
of the ionic Hubbard model The numerical re- 

sults are obtained with the DMFT (at zero temperature) 
and the determinant quantum Monte Carlo method. The 
existence of an intermediate metallic phase between the 
band and the Mott insulators seems confirmed by all the 
authors but the nature of this phase and of the metal- 
insulator transitions are still under debate. 
In this paper, we present results on the two-dimensional 
ionic Hubbard model obtained with a method based on 
the direct-space proposed by Suzuki and al.@, 0] and 
Hirsch and al.[?| Q. This quantum Monte Carlo method 
is presented in references [13, El ■ At fixed temperature, 
this method allows to generate some of the most repre- 
sentative occupation number basis states of the model. 
These states are used to compute average values of en- 
ergy, molar specific heat, occupancy, structure factor and 
rough static electric conductivity. 



II. IONIC HUBBARD MODEL 



The square lattice is a bipartite lattice with two sublat- 
tices A and B. c\ a and Cj, CT are the fermion creation 
and destruction operators at the lattice site i with spin 
c. i%i jC — cl a Ci jC is the number operator, t is the hop- 
ping term between nearest-neighbor sites, U denote the 
on-site Coulomb repulsion, A is the staggered potential 
between the A and B sublattices. 



III. SIMULATION PARAMETERS 

The square lattice contains 6x6 sites with periodic 
boundary conditions. Each elementary cell contains two 
sites A and two sites B. There are 18 spin up and 18 
spin down (half filling). The hopping parameter t is 
fixed at a value 1 except for the simulations of the 
atomic limit where t = 0. The model is decomposed 
in sub-systems which contain four sites. These sub- 
systems are grouped together in two sub-hamiltonians. 
The imaginary-time interval is divided into twenty slices. 
For each simulation five decreasing-increasing temper- 
ature cycles or decreasing-increasing interaction U cy- 
cles were programmed. There are one hundred points by 
curves. 



The Hamiltonian of the ionic Hubbard model can be 
written 



H = -t ( '•-.:'•/ ' + hc ) + U n ^ Ui t 



(i) 



iEB 



IV. THE ATOMIC LIMIT (t = 0) 

Simulations were performed for the simple case of the 
atomic limite where t — 0. For this value the Hamilto- 
nian is diagonal, so there is no problem due to the non- 
commutativity, in consequence the number of slides can 
be one, and there is no sign problem. Cycles of increas- 
ing and decreasing values of U were programmed at the 
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FIG. 1: (Color online). The energy at kT = 0.01 for different 
values of A. 
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FIG. 2: (Color online). The molar specific heat at kT = 0.01 
for different values of A. 



fixed temperature kT — 0.01 for different values of the 
staggered potentiel A. The conductivity is always zero, 
so the model is an insulator whatever the values of the 
interactions. The Figs. [TJ [5] and [3] display the results. As 
it is expected, one remarks that the electronic transition 
happens for U c = 2A without hysteresis phenomenon. 
There is one spin by site for U > 2A whereas only the 
sites of the sublattices B are occupied for U < 2 A. In 
this domain the energy of the model is E ~ 18 (U — 2A), 
while it is zero for U > 2A. There is not magnetic order. 

The Figs. HI [5] and M show the influence of the tem- 
perature on the energy curves, the molar specific heat 
curves and the occupancies curves versus interaction U. 
The occupancies curves for the different temperatures 
have very little error bars and cross almost exactly at 
U = 2A. This is in good agreement with the zero value 
of the specific heat for this value of U. At this point the 
site A occupancy is near 0.66 while the site B occupancy 
is about 1.32. In this special state 12 sites A are each 
occupied by one spin, 12 sites B are occupied equally 
by one spin and 6 sites B are occupied by two spins. 
The energy of this state is E = 6C7 — 12A = 0. Indeed, 
the three energy curves in Fig. 3] pass through the same 
point (U = 1,E = 0) so dE/dT = 0. 
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FIG. 3: (Color online). The site occupancies at kT = 0.01 for 
A = 0.5, A — 2 and A = 4. The solid lines corespond to sites 
A and the dashed lines correspond to sites B. 
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FIG. 4: (Color online). Energy in the atomic limit 
(t = 0, A = 0.5) for different temperatures. 



One remarks that all the specific heat curves of the 
Fig. [S] match exactly for the abscisse U/kT. The specific 
heat curves for four sizes of the model, at half filling, 
are shown in Fig. [7J One remarks that these curves are 
similar. 



V. MODEL WITH HOPPING INTERACTION 

(t = l) 

A. Simulations at T = constant 

The Figs. [SJ M EH and HU show the DC conductivity, 
the structure factor and the double occupancy curves for 
different values of the staggered potential A at kT — 
0.01. 

For A > 1 the model undergoes a transition between 
two insulator states. During the transition the system 
becomes conductor. In the low U insulator state, only the 
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FIG. 5: (Color online). Molar specific heat in the atomic limit 
(t — 0, A = 0.5) for different temperatures. 



FIG. 7: (Color online). Molar specific heat in the atomic limit 
(t — 0, A = 0.5) for four model sizes. 
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FIG. 6: (Color online). Site occupancies in the atomic limit 
(t — 0, A = 0.5) for different temperatures. Solid line (black) 
correspond to kT — 0.005, the dashed line (red) correspond 
to kT = 0.01, the dotted line (green) correspond to kT = 0.02 
and the dot-dash line (blue) correspond to kT — 0.05. 



B sites are occupied by two spins. It is a band insulator 
(BI). All the sites are occupied by one spin in the hight 
U insulator state. This last insulator state presents an 
antiferromagnetic structure, it is a Mott insulator (MI). 
This transition between the two insulator states with an 
intermediate metallic phase was already observed [IHEl 
with other simulation methods at T = and T ^ 0. Our 
results are in good agreement with those obtained within 
the other methods. 

Fig. [T^] shows the DC conductivity for A = 1 and A = 2 
for decreasing and increasing values of U. An hysteresis 
phenomenon appears for the transition from the metal to 
the Mott insulator while the the curves fit for the metal- 
band insulator transition. This behaviour is observed for 




FIG. 8: (Color online). Conductivity for low and hight values 
of U and different values of A. (t = l,kT = 0.01) . 




FIG. 9: (Color online). DC Conductivity for low values of U 
and different values of A. (t = 1, kT = 0.01) . 
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FIG. 10: (Color online). Structure factor for different values 
of A. (t = l,kT= 0.01) . 



FIG. 12: (Color online). DC conductivity for A = 1 and 
A = 2 for decreasing and increasing values of U at kT = 0.01. 
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FIG. 11: (Color online). Double occupancy for different val- 
ues of A. (t = l,kT = 0.01) . 
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FIG. 13: (Color online). Conductivity curves at night tem- 
peratures for A = 0.5 and different values of the coulombian 
repulsion U. 



all the values of A. One can deduce that the Ml-to- metal 
phase transition is a first order transition, while the BI- 
to-metal phase transition is continuous. This is in good 
agreement with the result of reference ■ 
One observes that the structure factor decreases for large 
A while the behaviour of the double occupancy is similar 
for all the values of A > 1.0. The caracteristcs of the BI 
phase ( null structure factor and double occupancy ps 0.5) 
are the same for all values of A whereas the insulator 
phase induce by increasing value of U is not a purely MI. 
Moreover, the conductivity of this phase is not null for 
A > 5. 



B. Simulations at U = constant 

The Figs. [T3] and [TJ] show the conductivity curves 
at low and hight temperatures for A = 0.5 and differ- 
ent values of U. One observes metallic behaviour for 



kT > 0.1. At low temperature, for decreasing temper- 
ature (Fig H"4"f , the system becomes insulator with be- 
haviour change for U ~ 1.75. This behaviour change can 
be observed equally on the double occupancy curves of 
Fig. [15] For U > 1.75 the metal-insulator transition 
occurs with hysteresis phenomenon (Fig. I16j) . One can 
deduce that this transition is a first order transition. On 
the contrary, for U < 1.75 the conductivity curves for 
increasing and decreasing temperature are similar. 



C. Phase diagram 

The conductivity curves of Figs. [8] and |9] can be used 
to drawn the phase diagram at the constant tempera- 
ture kT = 0.01. For large A one can consider that the 
metallic region shrinks to a single metallic point. For 
each value of A < 5.0, the coulombian interactions U c \ 
and U C 2 which correspond with the metal-insulator tran- 
sitions are determined at mid-height of the maximum 
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FIG. 14: (Color online). Conductivity curves at low tem- 
peratures for A = 0.5 and different values of the coulombian 
repulsion U. 
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FIG. 15: (Color online). Double occupancy curves at low tem- 
peratures for A = 0.5 and different values of the coulombian 
repulsion U. 
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FIG. 16: (Color online). Conductivity curves at low temper- 
atures for increasing and decreasing temperature for U = 2.0, 
U — 2.5 and U = 3.0 (A = 0.5). One observes an hysteresis 
phenomenon. 




FIG. 17: Phase diagram of the 2D IHM at kT = 0.01. 
The dot-dash line is the transition line at the atomic limit 
(U c = 2A). The symbol plus is the point at which behaviour 
changes. The dotted line correspond to the set of points for 
which double occupancy is 0.25. 



conductivity. The phase diagram of the ionic Hubbard 
model is shown in Fig. [T71 

The behaviour change observed on the double occu- 
pancy curves of Fig. [15] corresponds to the point 
(A = 0.5, U ~ 1.75) in the phase diagram. This suggests 
that a cross-over line exists in the metallic region. This 
line can correspond approximatly to the set of points for 
which the double occupancy is 0.25. For this double oc- 
cupancy value 9 sites A and 9 sites B are occupied by 
one spin and 9 sites B are occupied by two spins. 
One remarks, on the phase diagram, that the transition 
lines for t = and t = 1 are parallel for large A. The 
gap between these two lines is AU « It. 



VI. CONCLUSION 

Most results presented in references (H-Q are obtained 
for T — whereas, by principle, our simulation method 
works for not null temperature. However the results and 
the phase diagram are similar. Our results confirm the 
existence of a metallic region between Mott and band in- 
sulator phases. The natures of the metal-insulator tran- 
sitions are different. The MI- metal transition is discon- 
tinuous while the Bl-metal transistion is continuous like 
it is told in reference Q. The metallic phase shrinks to 
a line for large coulombian interaction U, but the BI-MI 
transition is not direct. Moreover, the insulator phase 
for U > U c is not a purely Mott insulator phase. Studies 
with increasing and decreasing temperature show that 
there is a behaviour change in the metallic region which 
divides it into two regions. These two regions correspond 
to the precursor phases of MI and BI phases. 
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